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In this paper we study the effect of a fluctuating gap in mono- and bilayer graphene, created by 
a symmetry-breaking random potential. We identify a continuous symmetry for the two-particle 
Green's function which is spontaneously broken in the average two-particle Green's function and 
leads to a massless fermion mode. Within a loop expansion it is shown that the massless mode is 
dominated on large scales by small loops. This result indicates diffusion of electrons. Although the 
diffusion mechanism is the same in mono- and in bilayer graphene, the amount of scattering is much 
stronger in the latter. Physical quantities at the neutrality point, such as the density of states, the 
diffusion coefficient and the conductivity, are determined by the one-particle scattering rate. All 
these quantities vanish at a critical value of the average symmetry-breaking potential, signaling a 
continuous transition to an insulating behavior. 

PACS numbers: 81.05.Uw,71.55.Ak,72.10.Bg,73.20.Jc 

I. INTRODUCTION 

Graphene is a single sheet of carbon atoms, where the latter form a honeycomb lattice. Graphene as 
well as a stack of two graphene sheets (i.e. a graphene bilayer) are semimetals with remarkably good 
conducting properties P, 0, These materials have been experimentally realized with external gates, 
which allow a continuous change of charge carriers. There exists a non-zero minimal conductivity at 
the charge neutralitypoint (NP). Its value is very robust and almost unaffected by disorder or thermal 
fluctuations [1, i, i, S] . 

Many technological applications of graphene require an electronic gap to construct switching devices. 
A first step in this direction has been achieved by recent experiments with hydrogenated graphene 0] 
and gated bilayer graphene [1, [l(| ■ These experiments take advantage of the fact that the breaking 
of a discrete symmetry of the lattice system opens a gap in the elctronic spectrum at the Fermi energy. 
A symmetry-breaking potential (SBP) is a staggered potential in the case of a monolayer, which breaks 
the sublattice symmetry of the honeycomb lattice, or a gate potential that distinguishes between the two 
layers in the case of bilayer graphene, where the latter breaks the symmetry between the layers. With this 
opportunity, one enters a new field, where one can switch between conducting and insulating regimes of a 
two-dimensional material, either by a chemical process (e.g. oxidation or hydrogenation) or by applying 
an external electric field 

The opening of a uniform gap destroys the metallic state immediately. Thus the conductivity of the 
material would drop from a finite value of order e 2 /h directly to zero. In a realistic system, however, the 
gap may not be uniform after turning on the SBP. This means that only locally the material becomes 
insulating, whereas in other regions of the sample it is still metallic. The situation can be compared 
with a classical random network of broken and unbroken bonds. The conductivity of such a network is 
nonzero as long as there is a percolating cluster of unbroken bonds. In such a system the transition from 
conducting to insulting behavior is presumably a second order percolation transition (l3j . 

Disorder in mono- and bilayer graphene has been the subject of a number of recent numerical studies 
[iH [l5| and analytic calculations [H, [l?], EH- The results can be summarized by the statement that 
chiral-symmetry preserving disorder provides delocalized states whereas a chiral-symmetry breaking scalar 
potential disorder leads to Anderson localization, even at the NP. The breaks the chiral symmetry but 
still allows for delocalized states at the NP [18|,[l!|. In contrast to chiral-symmetry preserving disorder, 
a random SBP reduces the minimal conductivity and can even lead to an insulating behavior. 

In this article an approach will be employed that eliminates a part of the complexity of the tight- 
binding model by focusing on continuous symmetries and corresponding spontaneous symmetry breaking. 
This allows us to identify a massless mode in the system with a randomly fluctuating SBP. Using a loop 
expansion we study the scaling behavior of the model and derive a diffusion propagator for the asymptotic 
behavior on large scales. Our result implies a relation between the average two-particle Green's function 
and the product two average one-particle Green's functions in self-consistent Born approximation. This 
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is similar to the solution of the Bethe-Sal pete r equation, a self-consistent equation for the average two- 
particle Green's function (Cooperon) [2(il l2ll I22I I23I l24j . In addition to the latter we are also able to 
control the scaling behavior of all higher order terms in the loop expansion. 

Our approach provides also information about the effect of symmetry-breaking terms. It turns out 
that the latter create a finite length scale Ldiff , such that diffusion breaks down for length scales L larger 
than Ldiff- Another reason for the breakdown of diffusion is a vanishing spontaneous symmetry breaking. 
This happens when the average value of the SBP exceeds a critical value. In this case there is no drop of 
the conductivity but a continuous decay to zero, depending on the fluctuations of the SBP. 

This article is organized as follows. In Sect. [TT] the model and functional- integral representation of 
the Green's functions are introduced. The symmetries of the model are discussed in Sect. Mil Then 
an effective functional integral is constructed for the average two-particle Green's function (Sect. IIV[) 
and a saddle-point approximation is employed (Sect. IIV A[) . The invariance of the saddle-point equation 
of Sect. IIV Al under a continuous symmetry transformation requires the integration over a saddle-point 
manifold. This is discussed in detail in Sect. [Vj which includes the loop expansion (Sect. IV Aj) . The 
results of the loop expansion and its consequences for the transport properties in graphene are discussed 
in Sect. IVII Finally, we conclude with a summary of our results in Sect. IVIII 

II. MODEL 

Quasiparticles in monolayer graphene (MLG) or in bilayer graphene (BLG) are described in tight- 
binding approximation by a nearest-neighbor hopping Hamiltonian 

H = t r _ T ic\c r i + V r c) r c r + h.c. , (1) 

<r,r'> r 

where c\ (c r ) are fermionic creation (annihilation) operators at lattice site r. The underlying lattice 
structure is either a honeycomb lattice (MLG) or two honeycomb lattices with Bernal stacking (BLG) 
[ll|,[l2j]. There is an intralayer hopping rate t and an interlayer hopping rate t± for BLG. V r is either a 
staggered potential (MLG) with V r = m on sublattice A and V r = —m on sublattice B, or it is a biased 
gate potential in BLG that is V r = m (V r — —m) on the upper (lower) graphene sheet. These potentials 
obviously break the sublattice symmetry of MLG and the symmetry between the two layers in BLG, 
respectivly. A staggered potential can be the result of chemical absorption of other atoms in MLG (e.g. 
oxygen or hydrogen Q). The potential in BLG has been realized as an external gate voltage, applied to 
the two layers of BLG [8|. A consequence of the symmetry breaking is the formation of a gap A g — m in 
both systems: The spectrum of MLG consists of two bands with dispersion 

E k = ±yjm? + (? k , (2) 

where 

t\ = < 2 [3 + 2cos/ci +4cos(fci/2)cos(\/3fc 2 /2)] (3) 
for lattice spacing a = 1. The spectrum of BLG consists of four bands [12] with two low-energy bands 



£ fe - M = ± V 4 + 4/2 + ™ 2 - + (*1 + 4 ™ 2 K (4) 

and two high-energy bands 



E+(m) =±\Jel + t\/2 + m 2 + + (t\ + 4m^ . (5) 

The spectrum of the low-energy bands has nodes for m = where (0) vanishes. These nodes are the 
same as those of a single layer. For small gating potential we can expand E^(m) under the square root 
near the nodes and get 



E fe H~±f 2 +E t -(0) 2 
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with the same gap as in MLG. 

The two bands in MLG and the two low-energy bands in BLG represent a spinor-1/2 wave function. 
This allows us to expand the corresponding Hamiltonian in terms of Pauli matrices <jj as 

H = h\o\ + h 2 a 2 + . (6) 



Near each nodes the coefficients hj are [23j 

hj^iV-i {MLG), hi = V? - V%, h 2 = 2ViV 2 (BLG), (7) 

where (Vi, V2) is the 2D gradient. 

Neither in MLG nor in BLG the potential is uniform. The reason in the case of MLG is that fluctuations 
appear in the coverage of the MLG by additional non-carbon atoms. In the case of BLG it is crucial 
that the graphene sheets are not planar but create ripples 0, [U [H| . As a result, electrons experience a 
varying potential V r along each graphene sheet, and m in the Hamiltonian of Eq. ([6]) is random variable 
in space. For BLG it is assumed that the gate voltage is adjusted at the NP such that in average m r is 
exactly antisymmetric with respect to the two layers: (mi),„ = -(r»2)m- 

At first glance, the Hamiltonian in Eq. ([1]) is a standard hopping Hamiltonian with random potential 
V r . This is a model frequently used to study the generic case of Anderson localization [2{|. The dispersion, 
however, is special in the case of grap hene due to the honeycomb lattice: at low energies it consists of 
two nodes (or valleys) K and K' [23l Wf\. It is assumed here that weak disorder scatters only at small 
momentum such that intervalley scattering, which requires large momentum at least near the NP, is not 
relevant and can be treated as a perturbation. Then each valley contributes separately to transport, and 
the contribution of the two valleys to the conductivity a is additive: a = uk + ctk'- This allows us to 
consider the low-energy Hamiltonian in Eqs. (J6j> , ([7]), even in the presence of randomness for each valley 
separately. Within this approximation the term m r is a random variable with mean value (m r ) m = m 
and variance ((m r — fh)(m r i — fh)) m = g5 r r > . The following transport calculations will be based entirely 
on the Hamiltonian of Eqs. I©,©. In particular, the average Hamiltonian (H) m can be diagonalized by 
Fourier transformation and is 

(H) m = fcicri + k 2 cr 2 + fha 3 



for MLG with eigenvalues Ek = zt^/m 2 + fc 2 . For BGL the average Hamiltonian is 

(H) m = (fc 2 ~ k 2 )a\ + 2k]k 2 a 2 + fhcr 3 



with eigenvalues — ±\/m 2 + fc 4 . 

Transport properties of the model can be calculated from the Kubo formula. Here we focus on interband 
scattering between states of energy ui/2 and —u>/2. This is related to the zitterbewegung [3(J, which is a 
major contribution to transport near the NP. The frequency-dependent conductivity then reads [l6j 

*o(o>) = ~^ 2 ((^-u/2\rt\^/ 2 )) m , (8) 
where \<&e) is the Fourier transform of the wave function under time evolution exp(— iHt)\ 

\^ B )= e ( - lE - e)t \^{t))dt= e {tE - e)t e- lHt dt\t>(0)) 



= -i(H -E- ae)- 1 ^^)) = —iG(—E - ie)|*(0)) (9) 

with the one-particle Green's function G(z) = (H + z) _1 . In other words, the conductivity is proportional 
to a matrix element of the position operator (fc — 1, 2) with respect to energy functions from the lower 
and the upper band. The matrix element on the right-hand side is identical with the two-particle Green's 
function: 

($- w /a|^|# w / 2 )=X) r * rr a[ G! ro(-w/2-ie)Go r (w/2 + ie)] . 

r 

With the identity H = —<j n H T <j n , where n = 1 for MLG and n = 2 for BLG (cf discussion in Sect. IIIip . 
the matrix element also reads 

(*-a,/2|r2|$ u/2 > = -^2rlTr 2 [a n G^ (uj/2 + ie)a n G 0r (i J j/2 + ie)] . (10) 
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A. Functional Integral 

The two-particle Green's function on the right-hand side of Eq. (fit))) can be expressed, before averaging, 
as a Gaussian functional integral with two independent Gaussian fields, a boson (complex) field Xrk and 
a fermion (Grassmann) field (k — 1,2) and their conjugate counterparts Xrk and [31(: 



- G I rrlijj ,(z)G rlr>klk (z) = J Vr'j'VrjXrkXr'k' exp(-So(*))X>[«TO • (11) 
Sq(z) is a quadratic form of the four-component field <j) r = [Xrii Xrii ^r2, ^2) 

S Q (z) = - j ^ <t> r ■ (H + z)r,r'0r' (imz > 0) , (12) 

r,r' 

where the extended Hamiltonian H = diag(H, H T ) of So acts in the boson and in the fermion sector 
separately. The use of the mixed field <f> r has the advantage that an extra normalization factor for the 
integral is avoided. The matrix element in Eq. (|10|) reads now 

(*-«/2|^|* w /2> =J2J2 r k [(*0i*rjXrfcX0fc) - (-l)"(*Oi*rfeXrjXOfc)o] 

= -J2J2 r k [(Xrk%jV0jX0k)a - {-l) n {Xrj^rk^OjXOk)o\ (13) 
jjik r 



with 



.exp(-S (z))V[V]V[x] 



III. SYMMETRIES 

Transport properties are controlled by the symmetry of the Hamiltonian and of the corresponding 
one-particle Green's function G(ie) = (H + ie)^ 1 . In the absence of sublattice-symmetry breaking (i.e. 
for m = 0), the Hamiltonian H — h\&\ + has a continuous chiral symmetry 

H —>■ e aa3 He°" 73 = H (14) 

with a continuous parameter a, since H anticommutes with 03. The term ma^ breaks the continuous 
chiral symmetry. However, the behavior under transposition hj = —hj for MLG and hj — hj for BLG 
provides a discrete symmetry: 

H -> -a n H T a n = H , (15) 

where n = 1 for MLG and n = 2 for BLG. This symmetry is broken for the one-particle Green's function 
G(ie) by the ie term. To see whether or not the symmetry is restored for e — > 0, the difference of G(ie) 
and the transformed Green's function —a n G T (ie)a„ must be evaluated: 

G(ie) + a n G T {ie)a n = G(ie) - G(-ie) . (16) 

For the diagonal elements this is the density of states at the NP p(E = 0) = po in the limit e — * . Thus 
the order parameter for spontaneous symmetry breaking is po- 

Eq. (|10[) indicates that transport properties are expressed by the two-particle Green's function 
G(ie)G(—ie). Each of the two Green's functions, G(ie) and G(—ie), can be considered as a random 
variable which are correlated due to the common random variable m r . Their distribution is defined by a 
joint distribution function P[G(ie),G(—ie)]. In terms of transport theory, both Green's functions must 
be included on equal footing. This is possible by introducing the extended Green's function 

A r x f G(ie) \ (H + it V 1 (M 
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In the present case one can use the symmetry transformation of H in Eq. (|15[) to write the extended 
Green's function as 

nUA-( a ° 0Wa o 0\(H + ie V'^o \ 
G(J£) -1 V -a o ){0 ia n ){ H T + ie) [o ia n J ' 

The extended Hamiltonian H — diag(H, H T ) is invariant under a global "rotation" 

H^e S He S =H, S = ( ,° (18) 

with continuous parameters a, a'. The invariance is a consequence of the fact that H anticommutes with 
S. The ie term of the Green's function also breaks this symmetry. For aa! — —tt 2 /4 the diagonal element 
of G — e s Ge is proportional to the density of states po- Thus, the continuous symmetry is spontaneously 
broken for e — > if po is nonzero. In this case there is a massless mode. 

As a symmetry-breaking parameter, e generates a characteristic response of the system with long- 
range correlations when it is varied for e ~ . This is reministent of a weak external magnetic field 
in a (classical) ferromagnet, where the response to a change of the magnetic field creates a power-law 
magnetic susceptibility near the critical point. Moreover, if e is chosen as a space-dependent field e r , we 
can vary it locally and obtain a space-dependent response in form of correlation functions of the Green's 
functions. This allows us to study complex correlation functions by taking local derivatives of the field 
e r . 

Returning to the quadratic form in the action So(z) of Eq. (fT2|) . we notice that after the "rotation" 

of the diag(H, H T ) with e s off-diagonal block matrices are generated. These matrices should have 
Grassmann elements in order to have a quadratic form that has pairs of complex and pairs of Grassmann 
variables. Therefore, the parameters a and a' must be Grassmann variables: a = if) and a' = ip. 



IV. AVERAGED MATRIX ELEMENTS 



As an example, we need to consider the averaged matrix element of r\ in Eq. (|13p . Averaging Eq. (jlip 
over the Gaussian distribution of v r means replacing exp(— So) by (exp(— So)) m on the right-hand side 
of the equation. The latter can be written again as an exponential function (exp(— So)) m = cxp(— Si), 
where the new function Si contains also quartic terms of the field tf>: 

Si = -i^tfir ■ {H + z) r y4> rl + 3 2j(^ r -73<M 2 ■ (19) 

r,r' r 

Then it is convenient to transform the integration variables (Hubbard-Stratonovich transformation [3l| ) 
as 

( XrXr Xrf r \ _^ A _ ( Qr ©r \ / 9n N 

where Q r , P r are symmetric 2x2 matrices and r , Q r are 2x2 matrices whose elements are independent 
Grassmann variables. Now the correlation functions in Eq. (|13p can be rewritten as correlation functions 
in the new field Q r . Then the matrix element reads 

((*- w /2|r2|* w/2 ))», - -^EE r " [(( Q ^)rM^)o,k 3 )2 - (-l) ,l ((ea 3 ) rjfc (8 ( 73)o jfc )2] (21) 

9 j^k r 

with 




...exp(-S 2 (z))OT2?[Q] 



and 

S 2 (z) = -Trg(Q^) + ln[detg[(i?) m + z - 2 73 Q]] . (22) 

r.r f ^ 



G 



Trg is the graded trace 

e b 

Tr is the conventional trace, and detg is the graded determinant |18j : 



T*g((i ? I )=Tr.\-TrB. 



A. Saddle-point approximation 

The integration in Eq. (|21| can be performed in saddle-point approximation. The saddle point is 
obtained as the solution of 8S2 = 0. Assuming a solution of the form 

V m s 

Qo = -1 ^73 - —70 , (24) 

we obtain the parameters 77, m s from the saddle-point equation 

Qo = g((H) m + z - 2 73 Q )- 1 73 . (25) 

A consequence of the symmetry discussed in Sect. IIIII is that for z = the saddle-point equation is 
invariant under the global symmetry transformation Qo — > U~ 1 QoU, where U — e s of Eq. psp . This 
transformation creates the saddle-point manifold 

Q'r = -i~7 S U 2 r ~ ^70 , (26) 

where U r is obtained from Eq. (|18[) by replacing the transformation parameters a (a') by space-dependent 
Grassmann variables ip r (ip r ), respectively. The form of Q' r , which is dictated by the symmetry, implies 
for the action S2 on the saddle-point manifold that (i) the quadratic term vanishes and (ii) the remaining 
term becomes 

S' = lndetg((iT) m + m s a 3 + z + i-qU 2 ) . (27) 
This action contains the symmetry breaking field z. The matrix element of Eq. pip becomes 

9 r 

with 

(...)s> = J ...e- s 'v[U] = J ...e- s 'v[i;} . (29) 

There is no extra factor from the invariant integration measure when we replace T>[U] by T>[ijj\ (cf 
Appendix [X| . 

B. Evaluation of the scattering rate 77 

The saddle-point approximation of the average one-particle Green's function means 

(G(z)) m = ((H + z)- 1 )™ » ((H) m + m s a 3 + z + i-q)- 1 = G (z + i V ) , (30) 

which is often called self-consistent Born approximation [20(. The ansatz for a uniform saddle-point 
solution in Eq. (|24p leads to a shift of z as z — > ir/' = irj + z with 

rf + iz = rj'gl (31) 
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and a shift of the average mass m — > m + m s with 

m s = -mgl/(l + gl) . (32) 

The integral I reads 



1 = 2 G ,nW)«ffc/(27r)7(i^) 



which is in the case of MLG 



1 i 

I / {rf + (m + to s ) 2 + k 2 y x kdk = — In 

7T ./n 27T 



A 2 



77' + (m + to s ) 2 



(33) 



and in the case of BLG 



1 



arctan ( A 2 / ^/ ?/ 2 + (m + to s ) 2 J 



I null (r ) ' 2 + (fh + m s ) 2 + k 4 )- 1 kdk= v ; = (34) 

n J° 2TrJri' 2 + (fh + m s ) 2 AJ -q' 2 + {fa + m s ) 2 

for A ~ 00. 

A nonzero solution rj for z — requires gl = 1 in Eq. (|3ip . such that m s = — to/2 from Eq. (|32[) . Since 
the integrals I are monotonically decreasing functions for large to, a real solution with gl = 1 exists only 
for |to| < m c . For both physical systems, MLG and BLG, the solutions read 

r/ 2 = (to 2 -to 2 )9(to 2 -to 2 )/4 , (35) 

where the model dependence enters only through the critical average parameter m c : 



2 A 



2Xe~ n / 9 (MLG) 



' . (36) 

. .9/2 (BLG) 



to c is much bigger for BGL, a result indicates that the effect of disorder is much stronger in BLG. This 
is also reflected by the scattering rate at to = which is 77 = m c /2. 

V. INTEGRATION OVER THE SADDLE-POINT MANIFOLD 

The integration weight exp(— S') of the functional integral in Eq. (|29|) reads according to Eq. ((27|) 

cxp(-S") = detg (h +ie + irp 2 ^ 1 (37) 

with the nonlinear field 

U 2 = e 2§ = 1 + 25 + 25 2 
and Hq = (H) m + ra s (j^. We notice that 

1 + s + s 2 = (1 -sy 1 , 

since S = for I > 3. This enables us to rewrite the integration weight as 

exp(S") = detg (h + ie-irj + 2ir/(l - S)- 1 ) = detg ((1 - S)(H +ie- irj) + 2ir)) detg(l - S)- 1 



detg 1 - S(H +ie- irj){H Q + ie + ir t )- 1 ) detg(l - S)- 1 , (38) 



where we have used that detg(iJ + if + iff) = 1. This result is remarkable because (i) S appears only 
linearly in the determinants and (ii) the matrix in the second determinant is diagonal: 

detg(l -§) = - 2&V> P ) • (39) 

r 

With the expression 

SGq := (H +ie- ir))(H + ie + ir])- 1 = 1 2ir](H + ie + ir/)" 1 = 1 - 2ir]G (i{e + r])) , 
we can write, using the definition of the graded determinant in Eq. (I23[) . 

Gxp(-S') = detg (l - SSGoj) 1 - 2^v) = det (l - ^£0,11^1^0,22) ~' JJ(1 - 2^ r ) . 

T T 

SGq depends on e, 77 and satisfies for n = 1 (MLG) or n = 2 (BLG) 

a n 5G ,n(e,r])a n = er„(l - 2irjG Q si(ie + ir)))a n = 1 + 2ir]G 0t2 2(-i( - iff) = <5G , 2 2(-e, -rf) ■ 
This implies for the integration weight 

exp(-S") = det (l - ^h^h+y 1 J|(l - 2t/v^v) (40) 

with h± = 5Go,22(ie, i??), whose Fourier components are 

h± = CT T 2ir]G ,22(±it ± «??) = c ± 2i?7cr n G'o,ii( = Re =F «?7))cn 



2i?7 

(77 + e) 2 + h\ + h\ 



a T : ^ 2 , rrro (=R(f? + e)o-o + (-l) n (/ii<7i - ft. 2 a 2 )) 



2rj(e + rj) 



(r) + e) 2 + h\ + h 2 2 



2ir 1 (-l) r 
(r) + e) 2 + ft- 2 + hj 



^o± — , ^ , l 2 , , 2 (-feui +/t 2 <r 2 ) ■ (41) 



Eq. (|40p is probably the most compact representation of exp(— S"), and a corresponding simple visu- 
alization is that the lattice has isolated sites (due to - 0V" T o) or closed random walks of h + and pairs 
(due to •tph_'iph + ). The functional integration in Eq. (f2"9")l can now be performed by expanding the de- 
terminant det (l — iph_4>h+) of Eq. (f4"01) in powers of the Grassmann variables ip r and Vv. A nonzero 
contribution to the integral requires that the entire lattice is covered with products VvVv- This is quite 
different from the corresponding functional integral with respect to complex fields, where already a single 
term of the expansion gives a nonzero contribution. Consequently, the expansion must be organized in 
a specific way to control the integration over the Grassmann variables. This can be done in terms of a 
loop expansion of the action S' which is discussed in the next section. 

A. Loop expansion 

Starting from the expression in Eq. (|4"0|) 

det (l — iph-iphy 1 = cxp(— In det (l — xph-iph + ^j) 
we can expand the exponent with trace terms of growing size as 

In det (1 - xjjh-iph+) = - ^2 jTr [(tph-iph+) 1 ] = ^ jTr [(h^iph+tp) 1 ] . (42) 

i>i i>i 
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FIG. 1: a) Elements of the loop expansion for the action S' and b) for the two-particle Green's function 
(G rr i(z)G r i r (— z)) m . The dot corresponds with a simple factor ip r ^pr from Eq. (|39jl . whereas the loops with I 
corners correspond with the expansion term of order I in Eq. (|42[) . Only an even number of corners can appear 
in the loop expansion a) and each site must be visited twice by line elements in b), except for the end points r 
and r' , which are visited once. 



The trace terms can be visualized as closed polygons (loops) on the lattice with alternating ip and ip at 
the corners (cf Fig. QJi), where each term is normalized by the number of corners of the loop I. Inserting 
this in the functional integral of Eq. , all the loops can contribute with the condition that they cover 
partially the lattice with products ip r ip r . There are many graphically equivalent coverages (but with 
different values), as can be seen in Fig. [TJa: a square can either be a product of four I = 2 contributions 
or just one I = 4 contribution. This equivalence raises the question for the contribution(s) to a given 
graph with highest weight in the functional integral. A way to study this is a scaling analysis, where we 
analyse the change of the loop-expansion terms under a change of length scales. For this purpose it is 
convenient to choose the Fourier representation 

Tr [(h-iph+i)) 1 ] = j ' ... J Tr 2 [h^ M ip kl ^ k2 h +M ip k2 ^ k3 ■ ■ ■ /i-,k a ,_ 1 V>fai_i-k 2 i^+,ibi^k2i-ki] d 2 k x ...d 2 k 2l . 

It should be noticed that there are only 2Z — 1 integrations that affect the field ip and its conjugate, 
namely fci — ki, k% — k$, k%i — k%, since the sum of these variables gives zero. The integration over the 
remaining 21 th variable affects only the h's. Using Aj = kj — fcj+i with k 2 i+i = k\ we get 



Tr [(h^h + ^) 1 ] = / C All ... Au^Ai^Aa ' • ' Va 2! _^a 2I <J(Ai + . . . + A 2l )d 2 A 1 ■ ■ ■ d 2 A 2l (43) 
with the coefficient 

Ca u ...,A2i = J Tr 2( h -Ai+---+A 2l +k 1 h + ,A 2 +---+A 2l +k 1 ■ ■ ■ h-,A 2l _ 1 +A 2l +k 1 h +: A 2l +k 1 )d 2 k 1 . (44) 

These integral expressions contribute with different weight to the loop expansion of exp(— S'), depending 
on the number of corners I. In order to analyse the weights we can use the fact that Aj as well as 
are integration variables in the functional integral. This enables us to rescale them as 

Aj-a^Aj, Va, -fl-^A, (45) 

and use the integration symbols as before the rescaling. Then the scaling behavior of the general loop- 
expansion term in Eq. (|43[) is 



Ca 1 ,...,a 2! V'a 1 ^a 2 • • • tpA 2 i_ 1 ipA 2l S(Ai + . . . + A 2l )d Ai • • • d A 2l 
s 2i(2+a) s -2 f c s a u ...,sA 21 ^aAa 2 ■ • ■ ^^A^Al + . . . + A 2l )d 2 A 1 ■ ■ ■ d 2 A 2l . (46) 
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Next, the contribution of C s Aj ... S A 2! to the prefactor must be determined. For I — 1 we have A 2 = — Ai 
such that 

Tr 2 (h_ M ij kl - k2 h + . k2 'ijj k2 _ kl )d 2 k 1 d 2 k2 = / tp Al / Tr 2 (h_ M h + _ Al+kl )d 2 k 1 ip_ Al d 2 A 1 



= J ^>A 1 C Al ^ Al d z A 1 . (47) 

This expression rescales as 

ip k C k ip_ k d 2 k -> s 2+2a / VfeC sfc ^_ fc d 2 fc , 



where C sk « Co + s 2 /c 2 Cq . Now we choose a = — 2 which gives a prefactor 1 for the s 2 £; 2 Cq term. 

In general, for s < 1 the rescaling of the wavevector in Eq. (|45[1 has the effect that the integration is 
shifted to larger values in Aj (i.e. to shorter scales in real space). This is compensated by a prefactor 
in front of the integral. A prefactor smaller than 1 means that the integral contributes more on larger 
values of kj than on smaller values. In other words, the corresponding loop contributes more to shorter 
length scales than to larger ones. Since we are interested in large-scale properties, terms with prefactors 
smaller than 1 are asymptotically irrelevant for this regime. The scaling of the coefficient 



a 



sAi sA, 



J Tr 2 {h^ tSAl ^ — \-sA 2l +k 1 h +:SA2 ^ hsA 2! +fci ■ ■ ■ h- tSA2l _ 1+sA2l+kl h +tSA2l+kl )d 2 ki , 



for / > 2 can be studied by rescaling h±. Then we have for each factor h± iSAj ^ hA 2! +Aii 

ft±,.A J +...A ai +fc 1 - h ±M + s(Aj + ■■■ + A 2l )h' ±M + o(s 2 ) . (48) 

such that 

21 21 I 

C.V -A . ■-(•, ■ >X C ^ A n + £ C iu~* II ^ + " 

jl=l ji,...,j ( =l n=l 

Here it is important to notice that each Aj becomes a gradient term in real space, whereas a constant 
term in Aj is diagonal in real space. Therefore, at least every second factor A^ (i.e., either A's with 
j = 1, 3, 21 — 1 or j = 2, 4, 21) must be present, since otherwise multiple factors of ip r or ip r at the 
same site r would appear which gives zero due to the fact that these are Grassmann variables. Thus the 
leading behavior of the right-hand side of Eq. (|46]) under scaling is 

For a — —2 this means that only terms with / < 2 are relevant for s ~ 0. Moreover, the I = 2 term 
vanishes, since there are two contributions that cancel each other. This can easily be seen in real-space 
representation: 

Tr(h-iph+'ij)h-iph + ip) = ^ Tr 2{h-^ 1 -r 2 h + ^ 2 -r 3 h-^ 3 -r i h + ^ r ^ ri )^ r2 : ij} rz 'tj} r ^ ri . 

n,...,r4 

The leading non- vanishing term is of order s 2 . In this case, according to the gradient expansion, every 
second term is diagonal and reads 

T^/l-^+^-ra^-^+^-rJVVi^raVVa^r! + Tr 2 (ft_, ri - r2 /l + ,o/l-,r 2 -ri fr+,o)VV 2 VV2^ViVVi ■ 

After renaming the summation indices and exchanging of the Grassmann factors in the first term we get 
= [-Tr2{h-fih +t r- L -v 2 h-fih + ,r 2 -r 1 ) + Tr 2 (h- t r 1 -r2h+,oh- t r 2 -r 1 h + ,o)] VV 2 VV 2 VViVVi ■ 



2i(2+ a ) s i-2 / c Ai ,...,a 2 ,Va 1 ^a 2 • • • ^^a^Ai + • • ■ + A 2l )d 2 A, ■ ■ ■ d 2 A 2l . 
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Now we use the fact that h± = KoOo ( K i a i + ^2^2) hi Eq. (|41[) and get from the sum of the two trace 
terms zero. This result implies that the loop expansion is asymptotically dominated by the term in Eq. 
(H7|) (i.e. the loop with two corners) which give the propagator 

£ e -*-(^)~_i_. (49) 



-2 + C q 

Here C q can be expanded in powers of q (cf. Appendix [B| as 

C g = 2-^4(e + lV)+ (q 3 ) 
9V 

with the diffusion coefficient 

D := -^--^ J Tr 2 [G a ^ k (e + ^Go^k^i-e - V )]d 2 k\ q=0 . (50) 

Thus the propagator in Eq. (j49|) describes diffusion on large scales. The e term corresponds with the 
symmetry breaking parameter. The latter does not need to be a scalar but can be any symmetry-breaking 
tensor in the Green's function, provided it allows us to write the two-particle Green's function in the 
form of Eq. 1(170. 

The matrix element of Eq. (|28[) reads with these expressions and the substitution e — ► ioj/2 



a n 1 1 

«* w /al*il*- u /a»m = -o3- 



g itu/2 + Dq 2 



rfD 

2 



9=0 £IW 



We can also use the definition of D in Eq. ([5TJ)) . together with Eq. (TIT))) , to write 



and 



£>=^-($? n ,| r ||$o. r; , ) (51) 



<(^/2l^l<f-./ 2 )) m = -r^^k^ -^) , (52) 



where \$%) is the wave function of Eq. @, when the Hamiltonian is replaced by the translational- 
invariant Hamiltonian Hq. Moreover, the integration in Eq. (f50|) gives for A ~ 00 (cf Appendix |C|) 



which implies 



D = r, os (a = 1 for MLG, a = 2 for BLG) (53) 

(4r/ 2 + m 2 )7r v ; y J 



8an' 2 

({^,2\r 2 k \$- u/2 )) m 2 ,,' 2 , ■ (54) 



VI. DISCUSSION 



All our results are obtained for the charge neutrality point E = 0, for mono- and for bilayer 
graphene. The main results are given in Eqs. (|49]) . |52|) . (fSTj) . and |53|) : Eq. (|49[) connects the average 
two-particle Green's function with the two-particle Green's function of the average Hamiltonian. A 
special consequence is Eq. (f52"|) . which describes a relation between a disorder-averaged matrix element 
and the corresponding matrix element of the pure system. Eq. (|5ip connects the matrix element with 
the diffusion coefficient. And finally, Eq. (|53[) connects the diffusion coefficient with the one-particle 
scattering rate 77. 
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density of states: The average density of states is proportional to the diagonal element of the average 
one-particle Green's function ((H + ie)~ 1 ) m . The latter can be evaluated in saddle-point approximation 
from Eq. (j2"5"|) as 

(G(ie)) m nG (ie + irj) , (55) 

where the parameters r\ (scattering rate) and m s are determined by the self-consistent (or saddle-point) 
conditions of Eqs. (|3Tj) , ([32]) . We then obtain po ps rj/2ng, where the scattering rate ij is a function of g 
and to, according to Eq. (135p . The density of states has a semicircular form with respect to to 

P ° ~ 2^ = 4^^"^ ~ m 2 e(m 2 - to 2 ) , (56) 
where the radius of the semicircle m c is given in Eq. (|36j) . 

diffusion: Scattering by the random gap term leads to diffusion, as explained in the loop expansion of 
Sect. IV Al The diffusion coefficient D in Eq. (|53[) depends only on rj' . This corresponds with the simple 
physical picture that diffusion decreases with an increasing scattering rate. Diffusion breaks down when 
the symmetry is broken by the parameter e. This implies a maximal diffusion length L^is = \J2D je. 
The scale L e ff indicates that any symmetry-breaking term creates a finite diffusion length which limits 
diffusion to systems of linear size L e g. This length scale can be very large due to the large diffusion 
coefficient D in MLG 



1 ge*/s 

In the case of BLG, however, it is much smaller because of the stronger scattering rate n = m c /2 of Eq. 

idiff ~ \ ■ 

V 7T6 



matrix element: The averaged matrix element ((& UJ /2\' r k\Q-uj/2))m is an indicator of Anderson localiza- 
tion, since it diverges if the localization length is infinite. According to Eq. (|5"3)) . the states \&± w /2) are 
delocalized at ui = 0. On the other hand, the states are localized for w/0 with a decreasing localization 
length as one goes away from the NP. Such a behavior was also found for bond disorder in analytic [l6| 
and in numerical studies (l5| . 

relation between averaged and non-averaged Green's functions: In general, the average two-particle Green's 
function can be expressed by the function C q through Eq. (|49| . C q in Eq. ()B2|) is a function of the 
Green's functions Go(±f/), where the random Hamiltonian H is replaced by the average Hamiltonian Hq. 
Since the average Hamiltonian is translational invariant, the function C q can be easily calculated. This 
relation between the average two-particle Green's function and the self-consistent two-particle Green's 
function 

J2e- iqr Tr 2 [(G r0 (-ie)G 0r (ie)) m ] « 

can be considered as a generalization of the self-consistent Born approximation of the one-particle Green's 
function in Eq. (|55|) . Like in the latter case, the averaging process leads to a change of energies e — » r/' 
(i.e. a replacement of the symmetry-breaking parameter by the scattering rate). A consequence for the 
matrix element is Eq. (|52p . which means a simple scaling relation between the average matrix element 
and the matrix element of the average translational-invariant Hamiltonian Hq. (The scale 77', however, 
is not free but fixed by the disorder average through the saddle-point equation (|25|) .) This provides 
an interesting and useful relation between averaged and non-averaged Green's functions. Moreover, in 
the relation of the matrix elements there is an extra prefactor —r)' 2 /(Lu/2) 2 . This is important for the 
transport properties, since it provides the derealization of states at lo = and it cancels the factor to 2 
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disorder g disorder g 

FIG. 2: Scattering rate r\ and diffusion coefficient D for m = in the case of monolayer graphene (full curves) 
and bilayer graphene (dashed curves) versus the variance g of the random symmetry-breaking potential. The 
diffussion coefficient of bilayer graphene is so small (D ~ 2/tt) such that it cannot be distinguished from the g 
axis. 



in the conductivity of Eq. (jHJ). The relation in Eq. (|5l2|) can also be understood as a factorization of 
the averaged matrix element into a product of a power law (i.e. ~ uj~ 2 ) and a smooth scaling function 

i 2 m n ,\ri\& iv ,). 



conductivity: The conductivity of Eq. (JSJ is calculated from the matrix element in Eq. (|54p and gives 

■nyvq 1 + tyl ) a 
It is remarkable that rf drops out for m = which gives a frequency-independent result 



2 

a e 
7T ft, 



A frequency-independent conductivity was also found for a random vector potential [la ] . In the absence 
of disorder a constant <r(ui) was found, with a different value though [32|,[33j]. The difference is due the 
fact that the expression in Eq. ([8]) is only a contribution due to interband scattering from the total Kubo 
formula (for details cf Ref. [la])- 

DC conductivity: For u> ~ the parameter rj' is replaced by the scattering rate 77 of Eq. (|36[) . The 
resulting DC conductivity reads 

7T(477 + m z ) /l 7T \ / n, 

Our knowledge of the diffusion coefficient D in Eq. ([55)1 and the density of states po in Eq. (|56p allows 
us to evaluate the DC conductivity alternatively through the Einstein relation: 

/ — 2 \ 2 
a / _ m \ , e 



«r( W ~0) K ^«^^l~je(^- ft ') fc . 

This agrees with Eq. ((58|) . except for a constant factor. 

It is important to notice that the conductivity at fh — does not depend on the variance g of the random 
SBP. This indicates that this quantity is robust against random fluctuations in graphene. In particular, 
we could have started from the action in Eq. (|19[) and treated the interaction term in perturbation theory 
in powers of g to obtain the same result. This idea was indeed employed in Ref. 17| and gave the same 
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value for the minimal conductivity. However, it is not possible to obtain a non-zero critical value m c in 
the case of MLG, since all orders of the expansion of m c in Eq. give zero. This is one of the reasons 
why we have not used the perturbation theory in g here but the loop expansion of Sect. IV Al 

VII. CONCLUSION 

The physics of the random gap model is characterized by a discrete symmetry of the Hamiltonian 
and a continuous symmetry of the two-particle Green's function. For the disorder-averaged two-particle 
Green's function the continuous symmetry is represented by a fermionic degree of freedom. Since the 
symmetry is spontaneously broken, the resulting massless fermion mode controls the properties on large 
scales. An effective action is derived for the massless fermion mode and a loop expansion is employed to 
extract the dominant large-scale contribution. It is found that the shortest loops are in control of the 
large scales, leading to diffusion. An explicitly broken symmetry generates a finite diffusion length Ldiff 
such that diffusion is possible only on length scales less than Ldiff- 

Although our models of mono- and bilayer graphene share the same type of symmetries and symmetry 
breaking, the quantitative properties are quite different, since scattering is much stronger in bilayer 
graphene (cf Fig. [2|). For instance, the diffusion coefficient D is very large for monolayer graphene, 
namely 

D cx ge n/9 

for average symmetry-breaking potential fh — 0, because the low density of states at the neutrality point 
does not provide much scattering. This means that transport in monolayer graphene is practically ballistic 
if disorder is not too strong. In the case of bilayer graphene, however, scattering is much stronger because 
of a large density of states at the neutrality point, leading to a constant diffusion coefficient D ~ 2/w for 
fh — 0. This also implies a large diffusion length scale Ldiff for monolayer graphene since Ldiff oc \f~D. 

All physical quantities of our discussion (i.e. the average density of states, the diffusion coefficient, 
and the matrix element of the position operator) depend on the model parameters only through the 
one-particle scattering rate r\. An exceptional case is the conductivity for vanishing average symmerty- 
breaking potential which is independent of the model parameters at all and has the value e 2 /irh for 
monolayer graphene and 2e 2 /nh for bilayer graphene (up to a factor 4 for spin and valley degeneracy). 
This implies a frequency-independent microwave conductivity. On the other hand, an increasing average 
symmetry-breaking potential fh reduces continuously the conductivity as well as the diffusion coefficient. 
The continuous behavior of the conductivity with respect to gap opening is similar to a recent experimental 
observation by Adam et al. (34[ • 
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APPENDIX A: INTEGRATION OVER THE NONLINEAR FIELD 

We consider the matrix expansion at fixed site r: 

Q = Qn + Qia + ■ ■ ■ 
where {Qij} is a basis for the matrix Q. In the integral 

I\ = / f(Q 11,11, Ql2,12) Q21,12, ■■■) 

Q1141 is replaced by the nonlinear term Qn,u + 1 + Q12, 12*921,12 which is created by the diagonal matrix 
elements of U 2 from the saddle-point manifold. This leads to the new integral 

I'2 = I f(Qll, 11 + 1 + *9l2, 12*921,12) *9l2, 12> *321, 12) •••) • 
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An expansion in terms of the Grassmann variable Qi2,i2Q2i,i2 gives 

1-2 = J /(Qll.ll + 1, Ql2,12, Q21 : 12, ■■■) + J Q Vl.YlQ 21, VI !\Q\\,\\ + 1, Ql2,12j ^21,12, ■■■) 

The second term vanishes at the Q 11,11 integration boundaries. Moreover, the shift in the first term by 
1 can be removed, since the integration of Qn 11 goes from — 00 to 00. This gives 



I2 — J /(Qll.ll, Ql2,12, G/21,12, ■••) — Ii 



APPENDIX B: DIFFUSION PROPAGATOR 

C q is defined in Eq. ((ii)) : 



C q = J Tr 2 {h^ :k h +:k ^q)d 2 k = J Tr 2 [a - 2i77G ,22,fc(«e + * 7 ?)][ cr o + 2iT]Go : 22,k~q(-ie - iil)]d 2 k 

{2 + 2ir)Tr 2 [Go,22,k- q (-ie - it]) - G ,22,fc(«e + irj)] + Arj 2 Tr 2 [G ,22,fc(«e + iri)Go,22,k- q (-ie - ii])]} d 2 k 



= 2+ J {2ir]Tr 2 [G Q ^2M-q(-ie - irj) ~ G ,22,k( ie + irj)} + 4r/ 2 Tr 2 [G ,22,fc(«e + ir))G ,22,k-g(-^ ~ "?)]} d 2 k 

( B1 ) 

since the k integral is normalized. The Green's function reads 

G ,22,fc(«e + irj) = - r — \ h 2 + h 2 t^( e + v) _ hl<Tl + h2<T ^ • 
Using the saddle-point equation (|25[) with 77' = r\ + e, we have 

?7 = ±i#Tr 2 [Go,22,rr(±«?/)] • 

This implies 

Tr 2 [G ,22,rr(-i??') ~ G ,22,rr («/)] = / 9 , 

such that 

C q = 2 - ^- + 4r, 2 / Tr 2 [G 0) 22,fc(ir/)G ,22,fc- g (-«/)] ^ • (B2) 
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The second term can be expanded in powers of q: 



C q = 2 - ^- + Ar, 2 I Tr 2 [Go,22,fe(^')G ,22, fe (- irf)] d 2 k 



fJ 



+ 2r] 2 qt— / Tr 2 [G , 2 2AW)Go^k- q {-W)} d 2 k\ q=0 + o{q 3 ) . (B3) 



£ 

■dq 2 

Since Go satisfies the following relations 

G (irf)G (-ir)') = (nf + ho)-\-W + ho)- 1 = (v' 2 + h^ 1 

and 

G (irf) - G (-irf) = (W + ho)- 1 - {-irf + ho)' 1 = -2w/(r/ 2 + h 2 )- 1 
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we obtain 



Tr 2 G (iv')Go(-iv') = ^Tr 2 [G a {i V ') - G (-i V ')} 



This allows us to write for the third term in Eq. (lB3j) 

V / Tr 2 [G ,22Aiv')G ,22A-W)]d 2 k 



2 3 

= VTr 2 [G ,22^ , )Go,22(-^')]rr = 2i\Tr 2 [G oa 2H) - Go,2 2 (-w/)]rr = 4^ 

V 9V 



This gives 



9V 



2 2 9V' d 2 



+ o(g 3 ) . 



Then the prefactor D of the q^. term reads 



D := - 2 ^ 2 I Tr 2 [G a ,22,k(W)G ,22,k-q(-iri')} d 2 k\ q=0 



APPENDIX C: EVALUATION OF THE MATRIX ELEMENT OF \$°± irt , 



The matrix element with respect to the average Hamiltonian (H) m of MLG gives 



<^K|$ u _ v )=4(^ + m74) 



dk 



o + fh 2 /A + k 2 ) 3 2n 2Tr(r)' 2 + m 2 /A) 



for A ~ 00, and of BLG 



(^l^l < f°V> = 16(^ 2 +™ 2 /4) / 

Jo 



k 3 dk 



1 



(r)' 2 + m 2 /4 + fc 4 ) 3 2-7T tt (r/' 2 + fh 2 / 4) 
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